/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include <math.h>
#include <stdlib.h>
#include <string.h>

static double
logprating (unsigned int        rating,
            double            (*alpha) (unsigned int),
            double              logbeta,
            unsigned int        n_labels) __attribute((unused));

static double
logprating (unsigned int        rating,
            double            (*alpha) (unsigned int), 
            double              logbeta,
            unsigned int        n_labels)
{
  /*
   * ok we want to compute the logarithm of
   * 
   * exp (exponents[rating])     exp (exponents[rating] - max)
   * ----------------------- = --------------------------------
   * sum_i exp (exponents[i])   sum_i exp (exponents[i] - max)
   *
   * the RHS is more numerically stable
   */

  double beta = exp (logbeta);
  unsigned int maxindex = 0;

  for (unsigned int i = 1; i < n_labels; ++i)
    {
      if (-alpha (i) * beta > -alpha (maxindex) * beta)
        {
          maxindex = i;
        }
    }

  double max = -alpha (maxindex) * beta;
  double denominator = 0;

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      denominator += exp (-alpha (i) * beta - max);
    }

  return -alpha (rating) * beta - max - log (denominator);
}

static void
dlogpratingdalpha (double*             dalpha,
                   unsigned int        rating,
                   double            (*alpha) (unsigned int),
                   double              logbeta,
                   unsigned int        n_labels) __attribute__((unused));

static void
dlogpratingdalpha (double*             dalpha,
                   unsigned int        rating,
                   double            (*alpha) (unsigned int),
                   double              logbeta,
                   unsigned int        n_labels) 
{
  double beta = exp (logbeta);
  unsigned int maxindex = 0;

  for (unsigned int i = 1; i < n_labels; ++i)
    {
      if (-alpha (i) * beta > -alpha (maxindex) * beta)
        {
          maxindex = i;
        }
    }

  double max = -alpha (maxindex) * beta;
  double denominator = 0;
  double ddenominator_dalpha[n_labels];

  memset (ddenominator_dalpha, 0, n_labels * sizeof (ddenominator_dalpha[0]));

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      double d = -alpha (i) * beta - max;
      double expd = exp (d);
      denominator += expd;
      ddenominator_dalpha[i] -= expd * beta;
      ddenominator_dalpha[maxindex] += expd * beta;
    }

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      dalpha[i] = -ddenominator_dalpha[i] / denominator;
    }

  dalpha[rating] -= beta;
  dalpha[maxindex] += beta;
}

static double
dlogpratingdlogbeta (unsigned int        rating,
                     double             (*alpha) (unsigned int),
                     double              logbeta,
                     unsigned int        n_labels) __attribute__((unused));

static double
dlogpratingdlogbeta (unsigned int        rating,
                     double             (*alpha) (unsigned int),
                     double              logbeta,
                     unsigned int        n_labels)
{
  double beta = exp (logbeta);
  unsigned int maxindex = 0;

  for (unsigned int i = 1; i < n_labels; ++i)
    {
      if (-alpha (i) * beta > -alpha (maxindex) * beta)
        {
          maxindex = i;
        }
    }

  double max = -alpha (maxindex) * beta;
  double denominator = 0;
  double ddenominator_dlogbeta = 0;

  for (unsigned int i = 0; i < n_labels; ++i)
    {
      double d = -alpha (i) * beta - max;
      double expd = exp (d);
      denominator += expd;
      ddenominator_dlogbeta += expd * d;
    }

  return -alpha (rating) * beta - max - ddenominator_dlogbeta / denominator;
}
